Wave functions for the Schwarschild black hole interior 

Daniel CartirQ 

Naval Academy Preparatory School, 197 Elliot Street, Newport, Rhode Island 02841 

Gaurav Kliann£0 

Physics Department, University of Massachusetts at Dartmouth, North Dartmouth, Massachusetts 02747 

(Dated: February 7, 2008) 

Using the Hamiltonian constraint derived by Ashtekar and Bojowald, we look for pre-classical 
wave functions in the Schwarzschild interior. In particular, when solving this difference equation by 
separation of variables, an inequality is obtained relating the Immirzi parameter 7 to the quantum 
ambiguity 5 appearing in the model. This bound is violated when we use a natural value for S based 
on loop quantum gravity together with a recent proposal for 7. We also present numerical solutions 
of the constraint. 
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Ph. I. INTRODUCTION 

< 

■ One of the more interesting class of solutions found in the theory of general relativity is the black hole. Not only 
does it fire the public imagination - as seen by a survey of popular books of both the fiction and non-fiction variety - 

' but it also provides an arena to study what happens when general relativity and quantum mechanics are important. 
^ ' The strong gravitational fields that curve space-time enough to prevent escape from the black hole interior eventually 

■ lead to a singularity, where the classical theory becomes meaningless. Like the Big Bang singularity, a complete 

■ understanding of what occurs at this location depends on uniting classical and quantum ideas. It is thought highly 
^vq , unlikely that there is a complete breakdown in the equations of gravity. Instead it is hypothesized that quantum effects 

' provide a limitation on the magnitude of space-time curvature, and perhaps even allow the possibility of discussing 
\^ , the region beyond the classical singularity. Thus the simplest form of black hole, the Schwarzschild space-time, is the 

■ focus of a variety of explorations into quantized gravity. _ 

I One of the more developed techniques is that of loop quantum gravity l]. Unfortunately, this theory still has 
Q^. some outstanding issues, so a symmetry reduced version of loop quantum gravity has been developed, known as loop 
^ ' quantum cosmology Q,y{. To some extent, this allows a testbed for different incarnations of the full theory, and leads 
(2J[), to the discovery of features that are robust under modification of the exact quantization method. Many models of 
1^ • cosmological interest have been studied using the wave functions of the full theory. The basic method is to start with 
I the kinematic Hilbert space of the full theory, then reducing down to those states obeying a particular symmetry 
' in order to quantize the Hamiltonian constraint. Research in cosmological singularities has been fruitful, and shown 
|h [ that there are no difficulties in resolving the singularity or indeed, in evolving the wave function through it Q . Now 
work has begun on the singularity occurring inside the spherically symmetric black hole (see also 0). This paper 
will show explicit solutions for the quantum constraint in the Schwarzschild interior and comment on the conditions 
necessary to ensure the wave function is smooth far away from the singularity, a region where quantum effects are 
not expected. This requirement captures the notion of pre-classicality Q, ensuring the wave function has the desired 
physical properties when there is not yet a physical inner product applicable to this situation. The expectation is 
that pre-classicality picks out wave functions in the eventual physical solution space. This is due to the fact that the 
Wheeler-De Witt equation for semi-classical states is an asymptotic limit of the quantum constraint equation. 

In Section^ we review previous work on the interior of the Schwarzschild black hole in the context of loop quantum 
cosmology. The interior portion of the space-time is chosen because it is of Kantowski-Sachs type. Since these metrics 
are spatially homogeneous, they are very similar to other models previously considered in loop quantum cosmology. 
The quantum operator corresponding to the Hamiltonian constraint is a partial difference equation acting on the 
eigenstates of the triad, and includes a quantum ambiguity S related to the fundamental length scale. Here, we point 
out the effect of using a self-adjoint constraint operator in the quantum theory. Because the wave functions will not 
have the same restrictions that are seen in other work (where a non-self-adjoint operator is used) the range of solutions 
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is correspondingly greater. This choice of constraint will eventually allow the study of the physical wave functions, via 
the use of group averaging. In SectionUm we use generating function techniques to solve explicitly for wave functions. 
To simplify this search, we use a separation of variables method to find solutions for the constraint equation. Solving 
for one of the two sequences is technically demanding, so we only consider its asymptotic limit in Section IIII CI 
deferring the exact results to the Appendix. In Section IIVI we report on numerical simulations of generic wave 
functions. Within the limits discussed in the paper, Gaussian wave packets can evolve under the quantum constraint, 
maintaining a classical trajectory until they get very close to the singularity. Imposing the physical condition that 
wave functions must be suitably smooth far away from the classical singularity leads to an inequality related the 
quantum ambiguity 6 mentioned above and the Immirzi parameter 7. Interestingly, this inequality is not satisfied for 
all the current proposals for 7 if a natural choice for S is made based on the full theory. These results are discussed 
further in Section Ivl 

II. THE SCHWARZSCHILD INTERIOR IN LOOP QUANTUM COSMOLOGY 

A model of the Schwarzschild interior in loop quantum cosmology has been carried out by Ashtekar and Bojowald 0, 
so we will briefly review its features here. Because this portion of the space-time is spatially homogeneous, the construct 
is similar to other models studied in LQC. The configuration space is coordinatized by the two independent triad 
components pb and pc, the line pb = gives the horizon of the black hole, while pc = is the singularity. The 
corresponding quantum operators for the triad act on their eigenstates by 

PfclM^T) = ^M7^pIm,t>. Pc|m,t) 7r^p|Ai,T), (1) 

where 7 is the Immirzi parameter, and £p is the Planck length. Thus the line n = locates the event horizon while 
r = is the classical singularity. There is a residual gauge symmetry on the phase space generated by the Gauss 
constraint of the model, which requires that a choice is made to fix this freedom. In previous work on loop quantum 
cosmology, a gauge choice such as p;, > was made. This still allows passage through the classical singularity, so its 
resolution can be studied. However, we can also consider sequences symmetric under ^ ^ —fi. This does not put an 
artificial restriction on the extent of the quantum wave function, but picks out those that satisfy the Gauss constraint 
of the model. We shall see in Section Hvl that this simplifies numerical simulations of the solution. 

The Hamiltonian constraint arising from loop quantum cosmology in the Schwarzschild interior acts on a triad 
eigenstate to give 

C^^^lfi, t) ^ {2-i''5Hl)-^[2{V^,+5,r - Vf,-s,rW + 2(5, T + 25) - l^i + 25, r - 2,5) - |/i - 25, r + 25) + \^i-25,T- 25)) 

+ {V^,r+5'V^,.r-5){\^i + ^S^T)-2{\ + 2-1^'5^)\^l,T)+\^l-A5,T))], 

where V^^^. is the eigenvalue of the volume operator, given by 

Notice that, because of the form of this constraint, the coefficients of the relation vanish at certain points. In particular, 
V^+s.T — y^i-s,T = when = 0, and V^,riA — T^,r-5 = when r = 0. This situation arises in previous work in loop 
quantum cosmology, such as the isotropic Q and the Bianchi I and IX models H, 0, . If we write the wave function 
solving the Hamiltonian constraint as a sum of eigenstates 

* '^^Sf,^r\^J.,T), 

then we find that 

2(V^-(5,T-25 — Vfi-3S^T-2S)Sp,-2S,T-25 + 2(Vfji,+zS^T+25 — Vfj,+5^T+2s)s^+2S,T+2S 

— '^{Vf_i-S,T+2S — ^^-3<5,r+2c5)s^-2c5,r+2<5 — 2(^^+35^.^-215 — V'^+i5,r-2i5)s^+2<5,r-2(5 
+ {Vfj,-4S,T+S — V^-4S.T-s)Sfj.-4S,T + {Vfi+4S,r+S — V^-|-45,r-(5)s^+4i5,r 

- (1 + 2j^5^){V^,r+S - V^,r-S)s^.r = 0. 

Thus so.r and s^^o never appear in any relation, since their coefficients always vanish in the recursion relation. 
As we will see below, the conditions for pre-classicality can be expressed in terms of relations between the sequence 
members s^.^ for low values of fj,, t. For a simpler example, we look at the recursion relation 

msm+i - 2Sm + msm-i = 0, m>l, 
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arising as a special case in the Bianchi I model Note that there is always a scaling freedom Sm — * cSm, so choosing 
the ratio Si/sq is the only free information we have in the initial values of the sequence. The sequence will oscillate 
without bound unless this ratio is fixed to be a certain non-zero value. Obviously this requires that both sq and si 
be non-zero. When we have vanishing coefficients in the recursion relation, as with the difference operator C^^\ a 
full range of smooth sequences is not available. This happens in the Bianchi I case, where sq — for exactly this 
reason. When we consider the self-adjoint version of the constraint, however, this situation does not arise and there 
are no vanishing coefhcients of the difference operator. So the boundary value sq is included in the partial difference 
equation coming from the Hamiltonian constraint and can be non-zero. 

This inclusion of the boundary values allows a much greater range of pre-classical sequences, as we will see later. 
However, in addition to this practical utility, the use of a self-adjoint version also helps in making contact with the 
physical Hilbert space of solutions. This is through the use of group averaging techniques 0, already used in a 
few applications to loop quantum cosmology p^ . Group averaging is a method to explicitly construct the physical 
solutions of a quantized and constrainted system. One picks out those functions solving the constraint by averaging 
over the one-parameter group generated by the constraint; this cannot be accomplished unless the constraint is 
self- adjoint. 

From this point on, we work with the symmetric operator Ti, — ^(C^^"> + C''^'^-'^). We write the constraint equation 

T-l\fi,T) = in terms of new parameters TO,n, where fi — 2m6 and r = 2nS. This results in the following recursion 
relation for all m, n: 



m 



\m+^\]V\n-l\ 



TO H — TO - 



m + — — m — — 



m H — — m 



^m— l,n— 1 
^m—l.n-\-l 
^m+l,n — 1 



(2) 



21 + \m\ 



1 



^m — 2,n 



(l + 2Y6')\m\{d\n+-\~J\ 



" + 51- 



n 

2 



Sm+2,n 



0. 



Because the original parameters /z, r can take any real value, there is no loss of generality by using the scaled variables 
TO, n. However, notice that the quantum ambiguity S and the Immirzi parameter 7 drop out almost everywhere when 
m,n are used, except for in the coefficient for Sm.n- This will be the only place where physical input will affect the 
solution. In addition, we note that the constraint greatly simplifies when \m\ > 3/2 a fact we will use in the next 
section. Now wc turn to the task of solving this difference equation. 



III. ANALYTIC SOLUTIONS 
A. Generating function techniques 

In general, it is difhcult to find solutions to multi-parameter recursion relations where the sequence is smooth for 
large values of the parameters. A generic feature are sign flips every time the sequence parameter is increased by one 
step. One effective means to study the behavior of these sequences is to employ a separation of variables technique, 
analogous to the procedure used for partial differential equations. The partial difference equation is then reduced to 
several one-parameter recursion relations, which are much easier to deal with. Examining the self-adjoint constraint 
for the Schwarzschild interior, we see that when we assume that m > 3/2 (i.e. consider the "bulk" of the space-time), 
the relation is separable in this manner. Denoting Sm,n = cumPm we get 

(to + l)a„i+2 - (1 + 27^(5^)toq:„i + (to - l)a„i-2 = A(q:„i+i - Q!,„_i) (3a) 
(/RrT^+/H)/3n+i-(yk^+/H)/3n-i = -A^yin+il-^ln-ilVn, (3b) 
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where A is a separation parameter. We can obtain solutions to these two recursion relations by using generating 
function methods which are reviewed later in this section. However, to use these techniques, the coefficients of 
the sequence values I3n must be polynomial. We factor out y/\n\ from the equation, and use a Taylor series expansion 
of the resulting coefficients, keeping only up to order 0(l/n) and considering only positive n for the moment^. This 
gives us 



Obviously this does not work when n = 0; when we look back at the original separated relation (|3b|) for /3„ with 
n = 0, we get the simple relation that (3i — (3-1. Finally, we can multiply our approximate relation Q for /3„ by 2n 
in order to give coefficients polynomial in n; because n > 0, this does not affect the resulting solutions. Thus, the two 
relations we will solve for separable sequences are the relation l|3a|) for am, and 

(1 + 4n)/3„+i + (1 - 4n)Pn-i - -A/3„ (5) 

for (3n ■ Because the recursion relation for the /3„ sequence is simpler ~ a second order relation as opposed to the fourth 
order difference equation for am - we will search for its solutions first. 

Solutions for a recursion relation in general will have undesirable physical properties. Generically these sequences 
will have oscillatory behavior - adjacent values will have the opposite sign, and the magnitude of the sequence values 
may increase without bound as the sequence parameter increases. At this time, there is limited work on finding a 
physical inner product in loop quantum cosmology, so there is no way to pick out such unphysical states. This leads to 
the notion of pre-classicality 0, where we put in place criteria that pick out wave functions with desirable properties. 
The expectation is that pre-classicality will pick out states that appear in the physical Hilbert space. This is based 
on the fact that the quantum Hamiltonian constraint is a discretized version of the Wheeler-De Witt equation ; in 
the limit of vanishing step size, solutions of the difference equation will go to those of the semi-classical differential 
equation. Thus, those sequences that have a smooth limit will match wave functions solving the Wheeler-De Witt 
equation. 

Because we are looking for sequences that represent the wave function of a space-time, which becomes classical for 
large volumes (i.e. far away from any singular points), we must have a way of restricting the wave functions to those 
that act semi-classically in the appropriate regime. With this in mind, the technique of using generating functions to 
solve difference equations has been developed in the context of loop quantum cosmology 10] ; for a review of these 
methods in a generic context, see Wilf Our goal will be to find a generating function B{y), whose Taylor series^ 
gives the sequence /3„, i.e. 

oo 

k=0 

Operations on the sequence /3„ can be mapped over to those on the function B{y). For example, suppose we want to 
find the generating function for the sequence /3„ = Pn+i, thus, we find that 

oo oo oo / oo \ 

k=0 fc=0 fc=l ^fc=0 ^ 

Thus the shifting operation /?„ Pn+i is equivalent to the operation B{y) [B{y)~ (3Q]/y on the generating function. 
Similarly, if we want a multiplication operator /3„ = nf3n, we get 

oo oc , oo 

k=0 k=0 ^ k=0 



^ Here we make a comment about the accuracy of using this approximate relation for I3„. When we look at the solution /3„ at large n, 
the errors between the order 1/n relation and the completely accurate equation will be very small, and limited to the regime where 
n ~ 0(1). Thus the difference between any solution we find for our relations here versus the relation for all m,n is noticeable only 
close to the singularity n = 0. Similar reasoning lets us ignore the difference between the full recursion relation and the "bulk" relation, 
where m > 2 allows us to simplify the absolute value signs for m. 

^ Notice that we are focussing solely on i nteg er values of the parameter fc, despite the fact that fi and r (and hence m,n) can take any 
real value. It has been shown elsewhere |l4j that using the sequence solution for integer values can be extended to all real numbers. 
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So multiplication /3„ njSn corresponds to using the Euler operator yd/dy acting on B{y). This is the reason we 
used the earlier approximation for the recursion relation; to get a differential equation for B{y) that is relatively easy 
to solve, we must stick to using operators that are polynomial in both the variable y and the Euler operator. 

Up to this point, the Taylor series is merely a formal device to get a solution to the generating function. However, 
we have also a physical requirement that the sequence is smooth for large values of the parameter, corresponding 
to classical space-times. This will correspond to limitations on the singularities of the function B{y), which can be 
seen in the following examples. Suppose we have a simple pole at ?/ = a, where < |a| < 1; the Taylor series of the 
monomial (1 — y/a)~^ alternates in sign if a < 0, and increases without bound when \a\ < 1, since 

Higher order poles makes the problem worse, since the magnitude of the coefficients of y" appearing in the Taylor 
series would increase accordingly. This gives an equivalence between singularities of the generating function B(y) and 
the asymptotic behavior of the sequence A singularity at y = 1, for example, shows this directly, since 

oc 

(1 - j/)B(y)|,=i = /3o + V(/3„+i - /3„)(1)"+^ = lim /3„. 

n=0 

When -B(j/) has a pole at y = 1 that is of order one or less, then the sequence /3„ approaches a finite value; if the pole has 
an order greater than one, the sequence is unbounded. Similarly, we can look for poles at y = — 1 to determine whether 
the sequence is oscillatory, hipping signs as the parameter n is incremented by one. If the generating function includes 
a monomial (1 + y)^^ having an order p close to one, oscillations may occur that eventually dampen out. In general, 
requiring that B{y) is finite whenever — 1 < j/ < will ensure the corresponding sequence /?„ has the appropriate 
semi-classical behavior. If we wish to have only bounded sequences, then we can add the further condition that B(y) 
is finite for < i/ < 1. Thus our criteria for the pre-classicality of a solution to the Hamiltonian constraint will be 
phrased in terms of the singularity structure of its associated generating function. 

We return to a simple example to illustrate what is happening, the previously mentioned relation 

ms,„+i - 2s„i + ms„i_i = 0, 

which arises from separation of the Hamiltonian constraint in Bianchi I LRS. The generating function for this relation 
is of the form 

, , flo — (2ao + a{)x — (4ao + 2ai) ln(l — x) _ C(x) 
^""^ " (TT^P " + 

where the function C{x) depends on two constants ao and a\ (these are the first two values of the sequence). Thus 
we have to worry about a second order pole at a: = —1, and we must pick a relation between oq and a\ to ensure the 
function F{x) is finite at this point. The reason we can do this is obvious when we Taylor expand C(x) around the 
point X = —\: 

C{x) =Co + C2(l+.x)^ + --- . 

Because C(x) has no term that is first order in (1 + then we can make cq = by the appropriate choice of constants 
ao and oi. Once this is done, the generating function Fix) is finite at a; = — 1, and the associated sequence will not 
have growing oscillatory behavior. If there was a non-zero term ci(l -I- x) in the expansion, the only solution would 
have been C(x) = 0, and the sequence would have zero for all its values. Once the pole at x = — 1 is taken care of, 
there is only the pole at a; = 1 due to the logarithm function. However, since (1 — x) ln(l — x) = at a; = 1, then 
limm^oo dm = implies the sequence asymptotically approaches zero. 



B. /3„ sequence 



In this manner we shall find a differential equation for the generating function, with the appropriate conditions on 
the order and location of any singularities. We make the shift i3(y) = /3o + yF{y) foi' convenience; from this we get 
the condition F{0) = Pi- In terms of F{y), the equation lO for /3„ corresponds to 



2,dF{y) 



42/(1 - y')^r^ + Fiy)il + Xy- ly') - 3/3oy - /3i = (7) 
dy 
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Solving this equation for F{y) gives 

F{y) = -iy-l/4(y_l)-3/4+A/8(^_^l)-3/4-A/8 



y 







(1 - j/)-3/4+A/8^^ _,_ j^^-3/4-A/8 



1 2 + A 2 - A 5 



3/3oy^ /5 2 + A 2- A 9 



^4' 8 ' 8 '4' 

5 ^'V4' 8 ' 8 '4'"""^"^ 

To ensure that -F'(O) = /3i, the integration constant obtained when integrating is set to zero. Fi is the AppeU 
hypergeometric function of two variables , defined as 

i'i(a,6i,f)2,c, zi,Z2) -2^2^ /„^ , -iLi ^1^2 

where (x)„ = r(a; + n)/r(a;) is the Pochhammcr symbol (or rising factorial), involving the gamma functions r(a;). 
Although an analytic continuation of this function can be defined for any values of the variables zi and Z2, as a 
series, Fi converges only when |zi| < 1 and |z2| < 1- AH of the singular behavior of the function B(y) generating the 
values of is contained in the function F{y). In this way, we can find all the properties of the sequences solving the 
relation ^ by studying the function F{y) solving a differential equation. Notice that when we solve for the function 
F{y), we find all values of the sequence simultaneously. Once the generating function is obtained, all values of the 
sequences can be read off the Taylor series expansion. This compares to using the recursion relation to find /3„, which 
would require first having Pn-i and Pn-2] the relation is inherently an evolution equation with a sense of "time" . The 
equation for Fly) is parametrized by the initial values (3o and pi, but it is not required we use these. Now that there 
is no restriction on Pq, as there would be in the non-self- adjoint case, we could use /3jv-i and Pn as our parameters, 
where iV is a large positive integer (see for more discussion). 

Now we locate the singularities of the function F{y); because of the polynomial appearing in front of the derivative 
in the differential equation, these are at y = ±1. Our only true degree of freedom is the ratio Pi/ Pq, since the sequence 
can be scaled by an arbitrary value, so we focus on keeping the generating function finite at y = — 1 to avoid growing 
oscillations. When evaluating F(— 1), there are the obvious poles due to the monomials in the function. For A < 6, 
there is a singularity at y = —1; similarly, there is another at y = 1 when A > —6. However, the AppeU functions also 
have their own singularities at y = ±1; to evaluate these, we use the fact [Tsj l that in the limit zi — > 1 

lim Fi(a; 61,62; c;zi,Z2) = M^(^^iji-^(1 - zi)^"'^-''^ (1 - ^2)"''^ (1 + 0(zi - 1)) 
21^1 1 (a)l (61) 

+ ^(^)^(^ ~ I'l ^FM b2;c - hi; Z2)(l + 0{zi - 1)) (9) 
i (c — a)l (c — 61) 

where r(n) is the gamma function and 2-^1 the Gaussian hypergeometric function. Recall that 2^1(07 h; c; x) is conver- 
gent only when the real part of (c — a — 6) > Il6| . Using the particular AppeU functions appearing in our generating 
function, and the exchange symmetry 

Fi{a;hi,b2;c; 21,2:2) = Fi (a; 62, 61; c; 22, 21), 

we find they are singular at y = 1 when A > 6, and at y = — 1 when A < —6. Thus, these functions have exactly the 
opposite behavior as the monomials. However, notice what happens as we vary A, for the y = — 1 pole in particular. 
When A < —6, the monomial (1 -I- y)~3/4-A/8 g^gg ^ero as y — > —1. As we can see from the behavior of the AppeU 
function as 21 1, and the exchange symmetry, the two AppeU functions diverge as (1 -I- y^/^+^Z^ . Thus, the two 
conspire so that the generating function as a whole remains finite at y = — 1. Because there is no singular behavior 
for A < —6, the associated sequence /?„ will not be oscillatory for any choice of values for /3o and Pi. However this 
does not occur when A > —6. In that case, the monomial now diverges; the AppeU functions are finite but non-zero, 
so there is no cancellation and F{~1) is divergent. Therefore the ratio Pi/ Pq cannot be arbitrary, but must be chosen 
by requiring 

/12 + A2-A.5. \ 3/3o /5.2-l-A 2-A 9. \ 
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This can be simplified by using the gamma functions of the second term in the expansion around zi = 1 after 
exchanging zi and Z2, giving 

/3i ( 6 \ 2Fi(t;i + i2+A;-i) 
/?o U + 8j2Fi(i;i + il + t;-l)' 

So, we find that the ratio Pi/po is completely free when A < —6, but must be fixed otherwise to avoid oscillations 
far from the singularity n = 0. We can expand the function B{y) — 0n + yF{y) in a Taylor series, to read off the 
values of the sequence /3„ for any real n, as done in previous work |l4|. However, in this case the function is not 
easily written in a compact form. Also, by examining the pole at y = 1 in a manner similar to the above, wc can 
see whether the resulting sequence /3„ is bounded or not. This gives us that the divergence of the Appell functions is 
cancelled out by the monomial going to zero when A > 2, so F{1) is finite in this regime, and the sequence goes to 
zero asymptotically. Otherwise, the sequence will increase without bound. Finally, we comment here that having a 
non-zero value /3o is crucial for obtaining pre-classical solutions for all values of the separation parameter A. If /3o = 
(as would be the case if a non-self-adjoint constraint were used), then the only sequence meeting the pre-classicality 
condition H1U|) would be the trivial one /3„ = 0. This would severely restrict the space of solutions, since physical wave 
functions would require A < — 6. 



C. Asymptotic limit of the a,n sequence 



When we treated the case of the /3„, the physical parameters did not enter anywhere in the analysis. Thus, the 
results there are independent of the Immirzi parameter 7 and the ambiguity S arising in the quantization. Because of 
the reparametrization we chose - going from the original triad eigenvalues fi and r over to the new variables m,n - 
the only place these appear is in the recursion relation for the separable sequence am in the combination ^S. As we 
will see, this will tie the existence of solutions with the proper semi-classical behavior to the values of 7 and S, and 
put a limit on their product. 

Because the study of the am sequence is more involved, we look here at the asymptotic behavior of the sequence, 
and reserve the full details to the Appendix. This will give us a roadmap to the results to be obtained for all values 
of the triad eigenvalue m using generating function techniques. By looking at the limiting cases of the sequences for 
large m, we will reproduce the division of behavior into unbounded and decaying sequences that occurs in the full 
solutions. To start with, we think of the recursion relation (|3a|) for am as an equation for a function a{ni). We use 
an expansion of a(m) in terms of a step size h, for example. 



^ a{m + 2h) = a(m) + 2 



da(m) 
dm 



h + 2 



d^a{ra) 
d'^m 



+ 0{h^ 



By doing the same for the rest of the relation, we obtain a differential equation for a(m) to various orders of h\ in the 
following, we go up to second order in h, and set h — \. This results in 



d'^a(m) ^^^da{m) . s. , . 

4m . \ ^ + (4 - 2A) ^ ' + (2 - K)a(m) = 0, 



dm? 



dm 



(11) 



where we define k = 1 + 27^(5^. Since the equation simplifies if we choose k — 2, let us solve it first for this case, 
giving 

a{m) = ao -t- [ai - ao]m^/^ + 0(m^/^~^), 

with the constants set using the first two values of the sequence am- Thus, already we can see what happens for large 
values of the triad parameter m when k = 2 ~ if A is positive, the sequence will increase without limit, while if A is 
negative, it will decay to zero. 

When K 7^ 2, the solution to the equation is 



a(m) 



,A/4 



A/4 



V2- 



C2^A/4 1 



(V2 



0(m 



A/4-l^ 



where J,y{x) and Y^{x) are Bessel functions of the first and second kind, respectively, and Ci, C2 are two constants of 
integration. Note that, when k > 2, the arguments of these functions will become imaginary, so the solution of the 
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equation will feature the modified Bessel functions I,y{x) and K„{x). We examine the large m behavior of k < 2 first. 
In this case, as a; — > oo, 

(^-) - --^/-^/^ COS - T - l) ' 

and similarly for Y^{x). Thus, we expect to find sequences that oscillate over a wavelength of 47r/-/2 — k, as it either 
decays (A < 2) or increases in amplitude (A > 2). Notice that our criterion for pre-classicality is to avoid the sequence 
changing signs as the parameter m increases by one step; a gentler sinusoidal oscillation is no grounds for discarding 
a solution since this often seen in quantum mechanics. On the other hand, when k > 2, the Bessel function I^{x) 
grows without bound for large x, since 

while the other function K^(x) exponentially decays. Since the exponential increases faster than any power of x, 
then the sequence will be unbounded for k > 2, regardless of the value of A. Thus, what we expect when we find 
complete solutions is that the am sequence will change its behavior as n crosses the critical value k = 2, from a slowly 
oscillating function to an unbounded one. 

IV. NUMERICAL SOLUTIONS 

Analyzing the separable solutions has enabled us to understand the general properties of pre-classical solutions 
of the Hamiltonian constraint. However, the structure of the self-adjoint constraint derived for this model is very 
beneficial in using numerical techniques. We now proceed to find solutions in this manner, using the parameter n (or 
r) as a time parameter. An arbitrary wave packet is chosen at some relatively large distance away from both the 
classical singularity and the horizon. Specifically we pick a profile Sm.N for a fixed large value N. This is done so 
that we can use the semi-classical approximation of the constraint, i.e. the Wheeler-de Witt equation given in 5], to 
find the derivative in n of the packet as it is evolved towards the singularity. 

Here we note several points that make the numerical calculation much easier. First, since we are using a self-adjoint 
constraint, there are no limitations on the values of the solution at the classical singularity. In previously considered 
cases, such as Bianchi class A models [T^ . those coefficients Sni,Ti2,n3 of ttic wave function corresponding to the zero 
volume basis elements drop out of the recursion relation, because of various factors of the volume eigenvalue appearing 
in the difference equation. For convenience in solving the equations in those models, the Sni,n2,n3 include a factor of 
the volume, which means that s„i.ri2,o — and similarly for the other boundaries. Meeting this requirement would 
mean the profile chosen far away from the singularity would have to be exactly right or else it would "miss" the 
correct boundary condition. This does not occur in the self-adjoint case - the relevant coefficients of the difference 
equation never vanish. Because there is no longer a restriction on the values at the singularity of the Schwarzschild 
interior, it is much easier to evolve arbitrary wave profiles. 

In a similar vein, the range of triad eigenvalues used to delimit the configuration space in the Schwarzschild model 
also results in an easier problem. When solving the constraints, there are some residual symmetries between the 
coefficients of the wave function The usual choice in previous loop quantum cosmology work has been to 

truncate the configuration space, and allowing only non-negative values of some of the coordinates. In the diagonal 
Bianchi class A models, for the eigenvalues of the triad operators, it is assumed that ni,n2 > (the third eigenvalue 
is unrestricted). Instead, here the wave functions are unaffected by the gauge transformation m — > — m, coming from 
the Gauss constraint 0- Putting this together, we can fix the boundaries of our numerical simulation to be m = ±M, 
for large M. This leads to the second simplification in numerically solving the equations. In the original work of 
Ashtekar and Bojowald, there is a discussion about the boundary conditions to be imposed on the wave functions. 
In particular, they make an argument for choosing Sm.n ^ as m — s- oo. The physical rationale for this choice is as 
follows. To make contact with classical general relativity, we want the ability to construct wave packets that represent 
a semi-classical wave function, peaked around the classical trajectories in the phase space. These are of the form 

Pbit) = pf^ ^t{2m~t), pS) = ±t\ (13) 

where m is the mass of the black hole and t an affine parameter, while pf'^ is a scaling factor that can be absorbed 
into the radial coordinate of the metric. We should not have solutions increasing monotonically for large pf,, since 
Pb{t) arcs back to zero as <: ^ 2m. Thus, as long as the maximum grid size M is large enough to avoid sizeable 
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errors in solving the difference equation, we can set s±M,n = 0. This gives all the information necessary to find a 
full solution for the Hamiltonian constraint. Here again, we see the problem of using k > 2; this boundary condition 
would be impossible to enforce with only unbounded sequences with A-dependent slopes. It can be done in the k < 2 
case, where all the sequences have the same asymptotic period of oscillation. 

When all these considerations are taken into account, one can obtain a numerical solution to the full Hamiltonian 
constraint, shown in Figure ^ In order to compare this to the classical case, we use the relations (jlBfl for the phase 
space trajectory, solve these for a relationship between the two momenta ph and pc, and compare this to the average 
value of the triad eigenvalues, /i and t (or equivalently, m and n), recalling the relations between the operator 
equivalents of pb andpc, and their eigenvalues /z and r. This comparison is done in Figure[21for a particular numerical 
wave function. As can be seen, the classical solution is a good approximation to the quantum wave function even very 
close to the classical singularity. Both the analytic and numerical sides have shown us a rich variety of wave functions 
that solve the Hamiltonian constraint for the Schwarzschild interior. From the form of the full constraint we can 
see that the wave functions will be symmetric on both sides of the classical singularity. This has implications for 
important issues in black hole physics, in particular, as seen in a recent proposal for a paradigm of information loss 
suggest by loop quantum cosmology results . 




FIG. 1: A numerical simulation of an m-symmetric Gaussian wave function as it approaches the singularity of the black hole. 
The black line in the middle represents the horizon and the classical singularity is to the upper right. 



V. DISCUSSION 

In the preceding work we used generating function techniques to analyze the behavior of solutions to the quantum 
Hamiltonian constraint in the Schwarzschild interior. This was done by using a separation of variables method, and 
studying the resulting one-parameter sequences. The sequence /?„ does not depend on any of the physical parameters 
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n 



FIG. 2: A comparison between the classical black hole and the quantum wave function. The green line indicates the classical 
trajectory given by the relations l|18|l . while the red line is the averaged middle point of a Gaussian wave packet numerically 
evolved using the quantum constraint 

7 and S of the theory, so has the same form regardless of what particular case we are dealing with. On the contrary, 
the combination k = 1 + j'^S'^ serves as a determinant for the asymptotic behavior (far away from the event horizon) 
of the am sequence. The goal is to use these separable solutions to assemble a semi-classical wave packet, which is 
peaked around the classical trajectory, so that for large m, n, the solution acts according to the equations of general 
relativity. Because these trajectories in phase space are entirely within a compact region, we expect that the wave 
function will act in a similar manner. With this in mind, whether or not the separable sequences go to zero as 
771,71 ^ oo is important. Since only am changes its properties with k, we focus on it for the remainder. 

If K > 2, then aside from two isolated solutions when A = ±Ao = ±\/k + 2 (see the Appendix), all solutions for am 
are unbounded sequences with an asymptotically exponential profile. This means that it is impossible to construct 
a generic wave packet to represent the black hole wave function for these values of k. There is no way to match a 
linear combination of these separable sequences to sum up to a constant value at large m because different values 
of A will give sequences with different slopes in the limit, as given by the asymptotic function H12(l for the k > 2 
sequence. On the other hand, when k < 2, there is a pre-classical solution for any value of the initial data and the 
separation constant A, so arbitrary wave functions are easily constructed. Hence we have found there is a relation 
between the values of the Immirzi parameter 7 and the ambiguity 6 arising in the quantization of the Schwarzschild 
interior; generic solutions will only exist when 

j6 < I/V2. 

From the full theory of loop quantum gravity, using the smallest eigenvalue of the area operator and comparing it to 
the Schwarzschild case, we find that S — 2^/3 5]. The inequality above then gives 

7ma. - ^ = 0.204124... (14) 

Thus, by the imposition of a boundary condition coming from the need to link up loop quantum cosmology with 
semi-classical general relativity, an explicit bound is placed on the Immirzi parameter. Although this is a tentative 
result - in that wave functions in the interior only are considered without reference to the outside - it is interesting 
to see that 7, a parameter somewhat analogous to the 6 parameter in non-Abclian field theories, can be bounded by 
physical considerations. 

The value of the Immirzi parameter is typically obtained by making contact with the Bekenstein-Hawking equation 
for the entropy of a black hole, that is Sbh = A/ A, where A is the area of the black hole event horizon. Since entropy 
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is obtained by the logarithm of the number of states contributing to a macroscopic state, the question is how to 
count these states. In loop quantum gravity, the area of a surface depends on the number of spin network edges that 
puncture the surface. Each edge is labelled by a spin representation {jk} G Z/2, for fc = 1, • • • , iV, giving a surface 
area of 

N 

A = 8njfp V\jk\{\jk\ + l). (15) 
fe=i 

Originally, the assumption was that for large surface areas, the lowest spin values would predominate the sum. 
However, this was shown to be incorrect by Meissner ■ In the limit of large area, counting all possible combinations 
of the spin labels give the equation 

oo 

k=l 

for the Immirzi parameter "/m- Solving this numerically gives 7m = 0.237532..., which violates the bound l|14|) 
obtained here. 

With this in mind, it is worthwhile to note some of the other proposals for the Immirzi parameter, where changes 
in the method of counting states is considered. For example, when it was assumed that the lowest order spin labels 
comprise the vast majority of entries in the sum, a correspondence between the area operator spectrum (|15f) and the 
classical quasi- normal oscillation modes of a black hole was observed [23|. The argument is the following. In the limit 
of large damping, the real part of the quasi-normal mode frequencies becomes lhqnm = ln3/87r. On the other hand, 
adding or subtracting a single puncture gives a change in area of 

AA = 87r7^2^|j™„|(|j™„| + l). (17) 

Equating these two results using Bohr's correspondence principle, one can show that "fQNM = ln3/27r\/2 = 
0.123637... and jVnm = 1- Although this is smaller than ^max, it requires that only integer spins are counted 
in the sum over states (due to the value oi jminY ■ There is the advantage of linking the quantum mechanical entropy 
calculation to a macroscopic quantity that could be measured soon. Thus, a mechanism for ensuring that all edges of 
the spin network puncturing the event horizon have equal values has been suggested 22j. 

If we take the Meissner value seriously, then one of the assumptions made to find the wave functions for the 
Schwarzschild interior is incorrect. Recall that we use the following ideas to obtain this bound: (1) pre-classicality to 
obtain sequences that are smooth far from the classical singularity; (2) the boundary condition Sm,n — > as to — > oo 
coming from the need to match with semi-classical physics; and (3) the choice oi 6 = 2^3, based on the smallest 
area eigenvalue of the full theory. Since assumption (1) is necessary to make contact with the Wheeler-De Witt 
equation, the continuous limit of the quantum Hamiltonian, and (2) is in a similar vein, it is unlikely that they are 
the culprits for this inconsistency. However, it is possible the choice of 5 is suspect, since it is based in part on the 
fiducial metric used to accomplish the quantization j5|, i23i | . In other words, the current method of quantization does 
not fix the physical value of the smallest area, since it depends on unmeasurable values - for example, the scale factor 
for an isotropic universe. An alternate method of quantization, where the area eigenvalues are determined only by 
physical quantities would result in a different value of 5 (in fact, it would be a function of the triad eigenvalues), so 
the inconsistency may not result in that case. Work on this is ongoing at the moment 12^; it remains to be seen 
whether the Immirzi parameter bound ^max will change enough to include the Meissner value 7m. 

The authors appreciate the helpful comments of Abhay Ashtekar, Martin Bojowald and the other members of the 
Institute for Gravitational Physics and Geometry in writing this manuscript. GK is grateful for research support from 
the University of Massachusetts at Dartmouth, as well as the Glaser Trust. 

APPENDIX: EXACT SOLUTIONS OF THE a™ RECURSION RELATION 

Since pre-classicality focuses on the behavior of the a,-n sequences for large eigenvalues of the triad, the analysis in 
Section inm reflects many important features. However, the sensitivity of the solutions to the choice of initial values 



^ Interestingly, if the full calculation is done using the Meissner relation I16i . using only integer spins, the numerical value 7 = 0.137727 
is obtained, again well within the bound. 
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also needs to be studied, by examining the sequences for all values of the parameters. Generic solutions of a difference 
equation will feature oscillations that increase without bound; these are the sequences we label as unphysical by using 
the notion of pre-classicality. To weed out solutions of this type, we have to find the sequence for all values of m by 
solving for its generating function A{x)] we will follow many of the steps we saw when working with the /?„ sequence. 
In particular, we start by defining a function F{x) associated to the full generating function A{x) by 

A{x) - ao - aix - a2x'^ - a^x^ 
F{x) = 2 . (A.l) 



Working from the recursion relation 13a|) . we arrive at the following differential equation for F( 



X 



dF(x) 

x(x'^ - Kx'^ + 1) + F{x)(5x'^ + \x'^ - Akx'^ ~ Xx + 3) + (ao + Xai - 2Ka2 - Aas) 

ox 

+ (2ai + Aq:2 - 3Ka3)x + (3a2 + Aa3)a;2 + 4q!3X^ = 0. (A. 2) 

The polynomial (x* — kx^ + 1) factors into monomials with roots at 

1 



xo = -(V« + 2 + Vk- 2), (A.3) 

as well ^, and — ^. The properties of these roots will play a role in determining whether pre-classical 

solutions are available for a particular k. Already, we can see that the properties of the polynomial appearing in front 
of the derivative of the generating function change at the value k = 2. We shall see that this is the point where range 
of pre-classical wave functions of the wave function changes as well. When we solve for the generating function F{x), 
we find that 

F{x) = -x-\l - xoa;)-i/2-V2Ao _ x/ xo)-^'^-'''^^'' {I + xox)-i/2+a/2Ao _^ x/xo)-^'^+'''^^" 

X j [(ao + Aai — 2Ka2 — Xa^) + (2ai + Aq!2 — 'iK,a^)x + (3a2 + Aa3)x^ + 4a3a;^] (A. 4) 

X (1 - Xoz)-^'^+^/^^°{l - z/xo)-l/2+A/2Ao (1 ^ ^^^^-l/2-A/2Ao ^ ^/a;o)-l/2-A/2Ao 



where the critical A is given by Aq = a;o + a^o ^ vk + 2. Obviously, the function F{x) must be the same regardless 
of which of the four roots we use. This can be seen by the invariance of Ao under the exchange xo ^ and why 

Ao ^ — Ao when xo —xo- The integral above is solvable in terms of Lauricella functions F'jy' of four variables [T5|. 
giving 

. X -l/2-A/2Ao / ^ X -l/2-A/2Ao / n -1/2+A/2Ao / ^ x -l/2+A/2Ao 

F{x) = -{l-xox\ V ~) [l + xox\ 



t rf.Fj^' + 3, + 4; ^,xox, -^,-xox ) . (A.5) 

^ fc + 3 V 2Ao 2Ao 2Ao 2Ao xo xo ) 

Here, the coefficients are 

Co = — A(q;3 — ai) — 2^02 + ao, c\ = — 3Ka3 + Aa2 + 2q!i, C2 = Aq!3 + 3q:2, Ci ~ 4a3 



and the integration constant has been chosen to be zero so that -F(O) gives the next value in the sequence, 014 = cq, 
after the four initial values; f'^^ is a four-variable extension of the Gaussian and Appell hypergeometric functions, 
defined by 

^(4). , , , , x_ (a)mi-|----|-m4(^l)mi • ■ ■ (64 



(a;fci, • ■ • ,^4;c;zi, . . . ,24) = X! "' X! 



'47m4 mi 



n n (c)mi-|---+m4Wl! • • • 7714! 

7ni— m4— V / I 4 

Analogous to the Appell function, the Lauricella function will converge as a series only when |zfc| < 1, for fc = 1, . . . , 4. 
This will become important in the k > 2 case. 

First we examine the situation where k = 2. This choice of k simplifies the generating function greatly, since 
Xo = Xq^ = 1 (giving Aq = 2), and the Lauricella functions are reduced to Appell hypergeometric functions: 

F{x) = (1 - :r)-i-V2(i + ^)-i+A/2 (^k + 3;l-^,l + ^;k + 4;x, -x^ . (A.6) 
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Now we have a case similar to the /3„ sequence, and we can use the expansion of the Appell function around the point 
a; = 1, given by @ and the exchange symmetry to see what happens at x = —1. Let us look first at what happens 
when X — —\. The monomial in front gives a pole of order at least one when A < 0, and is zero when A > 2. As we see 
from the first term in the Appell expansion (O, the hypergeometric function has singular behavior like (1 + a;)~^/^, 
giving a finite value for that term when A < 0. Up to this point, we have the same results with the /?„ sequence - 
any divergence in cither the Appell function or the monomial is balanced by the reciprocal in the other. However, 
here we have a new wrinkle in the second term. Previously, with /?„, the Gaussian hypergeometric function in the 
coefficient was finite, so we did not worry about it. This is not the case here, since we have for each term in the sum 
+ 3, 1 — A/2, k+'i — A/2, —x), which is divergent^ for all values of A This introduces the need for a relation 
between the to make the second term in the expansion around x = —1 finite, regardless of A. After this, we need 
a second relation to insure F{—1) is finite, just as in the /3„ case; it has taken two restrictions on the initial data to 
remove any oscillatory behavior in the sequence. When we turn to what happens at x = 1, we find that again, we have 
a divergent Gaussian hypergeometric function as we approach for all A. Thus our remaining freedom in the initial 
values is taken up ensuring the Appell function is not divergent. Once this is done, the first term in the expansion of 
the Appell function comes into play, with its singular behavior as (1 — a;)'*'/^. Then when A < 0, the Appell function is 
divergent, but together with the monomial, F{x) ~ (1 — x)~^ near x = 1 and the sequence is bounded. When A > 0, 
the monomial is singular but the hypergeometric function is non-zero (because of the second term in now finite 
by selection of the c^), so F{x) ~ (1 — x)~^~^/'^ , giving an unbounded sequence. 

Next we look at k 7^ 2. When k > 2, the root xq will be a positive real number greater than one, while it is a 
complex number with unit modulus if k < 2. In both cases, we can use the formula (|A.4|I given above, where there are 
potential singularities at ±a;o and ±a;(7^- Analogous to their simpler cousins the Appell functions, the parameters in 
the Lauricella function show there will be a pair of poles for a given value of A, either at xq and x^^ , or else —Xq and 
its reciprocal. For example, the hypergeometric function will have a pole at —x^^ (which gives unbound oscillations) 
of order 1/2 — A/2Ao, i.e. when A > Aq. Similarly, there is a pole at x = Xq^ when A < Aq. We immediately run into a 
problem with the roots ±xo, however, because the Lauricella functions as a series expansion is convergent only when 
the variables \zk\ < 1. Thus when k > 2, we have functions in the sum of the form 

~ ^ ^0 ~ ^ ^" + ^ ^0 + ^ ./, I j. 1 -I 2 

\ 

in the sum for F{—xo). A similar result is obtained at a; = xq- We saw when discussing the behavior of series, in 
particular the series © for a simple pole at y = a that there may be oscillatory behavior for a finite order pole in the 
interval (— cxd, —1), but it will eventually decay in amplitude. Here, the Lauricella functions will diverge faster than 
any finite power of x due to the two arguments that are outside the unit box \zk\ < 1. We can use the expansion 

(a) (6l)mi(^2)m2(^3) 



F^\a;bi,...,bi;c;zi,...,Z4:) = X! ' ' ' X! 



n n {c)mi+rn2+m3milm,2\m3\ 

mi— m3— ^ ' 

/ 3 3 

X z"'^z'^^z!!^'''^2Fi ( a + X! ^4; c + ^ nik; Z4 
^ fe=o fc=o 

to analyze what happens at these divergent points; the choice of which of the four variables to expand around is 
obviously symmetric. Note that the Gaussian hypergeometric function in the sum will have the same convergence 
properties, regardless of the values of ruk, based on the real part of (c — a — 64). The result of this is that we must 
ensure the generating function has a finite order pole at a; = —xq to avoid oscillatory sequences coming from the 
divergence of the Lauricella functions. This requires two relations on the initial data, since there are two of the four 
variables with a magnitude greater than one. The final relation comes from requiring the generating function is finite 
at a; = —Xq^. Once we have gotten rid of the divergent pieces, at this value of x we either have a pole in the monomial 
for A > — Aq, or else in the Lauricella function when A < — Aq. This is a situation similar to the /3„ sequence and 
the K — 2 with a similar resolution. Specifically, F{—X()) is finite if A > — Aq, but requires a extra condition on the 
initial data if A < — Aq (because the monomial is divergent but the Lauricella function is finite). It is important to 
note that there is not enough freedom to completely cancel out the unbounded rise of the sequence; evaluating F{x) 



It is important to realize that we are speaking here about divergence as a series. The hypergeometric function can be analytically 
continued so that a finite value is obtained at a; = — 1; this is what is done in mathematical software such as Maple and Mathematica. 
However, graphing the coefficients of series definition of the function will show that it has exactly the unbounded oscillatory properties 
that we want to eliminate. 
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at X = Xq gives the same types of divergences in the series, and one cannot impose enough conditions so their poles 
in the generating function are of finite order. 

However, there are two particular cases for n > 2 where it is possible to get a bounded solution - namely, choosing 
A = ±Ao - which are not obvious from the asymptotic analysis. If we choose A = Aq as an example, the generating 
function simplifies to become 



3 

k=0 



Ck Fi{k + 3, 1, 1, fc + 4, -x/xq, -xqx) 
(1 - a;o2;)(l - a;/a:^o) 



The Appell function will have poles at x — —xq and x = —Xq^, as can be seen from the discussion of these functions 
in the last section. Since the Gaussian hypergeometric function in @ is conditionally convergent at both of these 
values, there will only be two relations fixing the coefficients Ck (this can be seen more easily when the Appell function 
is written in terms of logarithms and polynomials). We do not have the A-dependent poles anymore, so this holds for 
any values of k > 2. For the remaining freedom of the coefficients Ck, we can choose to fix it so the generating function 



F{x) is finite at x 



to avoid unbounded solutions. This will give us a sequence that tends asymptotically to 



zero; the situation is similar in the case A = — Aq. These are isolated cases, however, and are not helpful in assembling 
a wave packet composed of a range of A. Here we have shown it is impossible to construct generic solutions for k > 2 
if one desires the wave function to have the proper semi-classical behavior far from the classical singularity. 

When K < 2, we do not have the same difficulty, since |a;o| = 1. So already we have more freedom in choosing 
initial values for the sequence. It turns out that the initial data are completely free, because of the complex nature 
of the roots. Specifically, we have xq = exp{i9), where 9 = cos~^{^/k + 2/2). For our range oi 1 < k < 2, 9 covers the 
interval < 6* < 7r/6, so none of the roots xq, —Xq, Xq^ and —Xq^ are identical. We can see what is happening in the 
sequence by looking at the Taylor series of the product of two monomials, whose roots are complex conjugates. This 
gives 



oo oo 



[(l-e^M(l-rt)f = EE 



j=0 k=Q 



3 J \k 



E E 

;=o '-m=0 

oo , [1/2] 

2E{E 

1=0 m=0 



P\ P \ .(/-2™)6 

m) \ l — ml 



P\l P 
ml \l — m 



cos[(/ - 2m)9)] }x 



with [1/2] the greatest integer less than or equal to 1/2. Obviously the binomials can cause the amplitude of the 
sequence to grow for the right range of p. Yet the fact we have cosine functions whose sign can change in the 
coefficients of cc' means that we get slower oscillations in the sequence, unlike the alternating sign changes of sequences 
not considered physical. Since this is a generic statement, it will hold regardless of the choices of initial values, so 
pre-classical sequences can be found for any choice. Two particular examples are shown in Figures 13 and ^ 
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